// STL
#include <iostream>                            // std::cerr, ::endl
#include <tuple>                               // std::tie()
#include <utility>                             // std::pair<>, std::make_pair()
#include <vector>                              // std::vector<>

// stats++
#include "statsxx/postprocess/ROC.hpp"         // ROC_pt, ROC_CostMatrix, cutoff_CB_well_behaved()


//
// DESC: Finds the optimal cutoff point on an ROC curve, using cost--benefit analysis.
//
// INPUT:
//
//     std::vector<ROC_pt> ROC     :
//     ROC_CostMatrix      cm      :
//
// OUTPUT:
//
//     bool                success : whether an optimal cutoff was found
//     double              cutoff  : optimal cutoff
//
// NOTE: A good discussion of cost/benefit analysis and ROC curves can be found in:
//
//     C. E. Metz ``Basic Principles of ROC Analysis'' (1978)
//
// NOTE: The operating point when the cost matrix has a slope of 1 is identical to optimizing the sum of sensitivity and specificity.
//
inline std::pair<
                 bool,
                 double
                 > cutoff_CB(
                             const std::vector<ROC_pt> &ROC,
                             // -----
                             const ROC_CostMatrix      &cm
                             )
{
    bool   success;
    double cutoff;

    // CALCULATE THE SLOPE AT THE BEST OPERATING POINT
    // TODO: NOTE: Maybe this goes to a (type of) utility subroutine.

    // NOTE: Return mathematical thresholds (0 or 1), rather than just above or below end-point scores of ROC.
    // NOTE: ... In these cases, one (likely) need not make predictions anyway.

    // for a prevalence P of 100% ...
    if( ((1. - cm.P) == 0.) || ((cm.FP - cm.TN) == 0.) )
    {
        return std::make_pair(
                              true,
                              0.
                              );
    }
    // ... else for a prevalence P of 0% ...
    else if( (cm.P == 0.) || ((cm.FN - cm.TP) == 0.) )
    {
        return std::make_pair(
                              true,
                              1.
                              );
    }

    double m = ((1. - cm.P)/cm.P)*( (cm.FP - cm.TN)/(cm.FN - cm.TP) );

    // CALCULATE SLOPE AT EACH POINT ON ROC
    // TODO: NOTE: Maybe this (also) goes to a (type of) utility subroutine.
    std::vector<double> slopes;

    for( auto i = 1; i < (ROC.size()-1); ++i )
    {
        double _slope = (ROC[i+1].TPR - ROC[i-1].TPR)/(ROC[i+1].FPR - ROC[i-1].FPR);

        slopes.push_back(_slope);
    }

    // NOTE: Slopes are added at the front and back similar to a natural spline (zero second derivative).
    slopes.insert(slopes.begin(), slopes.front());
    slopes.push_back(slopes.back());

    // DETERMINE IF ROC CURVE IS WELL BEHAVED (MONOTONICALLY-DECREASING SLOPE)
    // TODO: NOTE: Maybe this also goes to a (type of) utility subroutine.
    // TODO: NOTE: ... In this case, use the appropriate terminology --- whether the curve is "proper".
    bool well_behaved = true;

    // NOTE: Ignore the end points, where zero second derivatives have been assumed (so slopes are the same).
    for( auto i = 2; i < (slopes.size()-1); ++i )
    {
        if( slopes[i] >= slopes[i-1] )
        {
            well_behaved = false;

            break;
        }
    }

    // MATCH SLOPES
    if( well_behaved )
    {
        std::tie(
                 success,
                 cutoff
                 ) = cutoff_CB_well_behaved(
                                            ROC,
                                            slopes,
                                            // -----
                                            m
                                            );
    }
    else
    {
        std::cerr << "error in cutoff_CB(): code for non-well-behaved ROC curves is commented out" << std::endl;

        success = false;
    }

    // RETURN
    return std::make_pair(
                          success,
                          cutoff
                          );

/*
    // for a non well behaved curve, find the point on the ROC curve with the shortest line from it that intersects perpendicularly with a line with slope m passing through (0,1)
    else
    {
        // place a line with slope m that intercepts (0,1): y = m*x + 1

        // the line perpendicular to this has a slope of (-1/m)
        double minv = -1.0/m;

        decltype(ROC.size()) optidx;

        double min_d = std::numeric_limits<double>::max();

        for( decltype(ROC.size()) i = 0; i < ROC.size(); ++i )
        {
            // find the intercept of the equation with slope (-1/m) that passes through this ROC point: y = minv*x + b
            double b = ROC[i].TPR - minv*ROC[i].FPR;

            // now equate the two equations to find the x position of intersection: (m*x + 1 = minv*x + b); [(m - minv)*x = b - 1]; [x = (b - 1)/(m - minv)] ...
            double x = (b - 1.0)/(m - minv);
            // ... and use the first equation to get y
            double y = m*x + 1.0;

            double d = (x - ROC[i].FPR)*(x - ROC[i].FPR) + (y - ROC[i].TPR)*(y - ROC[i].TPR);

            if( d < min_d )
            {
                optidx = i;

                min_d = d;
            }
        }

        return std::make_pair( true, ROC[optidx].score );
    }
*/
}
